# Author: Mark Richardson
# Purpose: Format replication data for partisan disagreement

# Load packages
library(dplyr)
library(rstan)
library(stringr)

# Load data
load("data/01_disagreement_model_fitted.RData")

#### Format agency and mission disagreement for export ####

# Rescale to be N(0,1)

# Create variables for export from posterior draws

x <- as.matrix(pty_dis_fit, pars = c("mu", "theta"))

x_std <- (x - mean(x)) / sd(x)

mean(x_std) # Accuracy check
sd(x_std) # Accuracy check

x_out <- tibble(mean = apply(x_std, MARGIN = 2, FUN = mean),
                sd = apply(x_std, MARGIN = 2, FUN = sd),
                lb_95 = apply(x_std, MARGIN = 2, FUN = quantile, probs = 0.025),
                ub_95 = apply(x_std, MARGIN = 2, FUN = quantile, probs = 0.975),
                par = colnames(x_std))

mu_std <- x_out %>%
  filter(str_detect(par, "mu")) %>%
  mutate(mission_index = as.integer(str_remove(str_sub(par, start = 4), "]"))) %>%
  full_join(., mission_index, by = "mission_index")

theta_std <- x_out %>%
  filter(str_detect(par, "theta")) %>%
  mutate(office_index = as.integer(str_remove(str_sub(par, start = 7), "]"))) %>%
  full_join(., office_index, by = "office_index")

rm(list = ls()[!str_detect(ls(), "mu_std|theta_std|sfgs|x_std")])

# Apply population and respondent thresholds

thresh <- sfgs %>%
  mutate(pty_dis_complete = if_else( (!is.na(pty_dis_freq_int) & !is.na(pty_dis_strg_int)), 1, 0)) %>% # Indicator for complete set of responses
  group_by(office) %>%
  summarize(pop = n(),
            pty_dis_n = sum(pty_dis_complete)) %>%
  ungroup()

n_min <- 4

pop_min <- 20

theta_std <- left_join(theta_std, thresh, by = "office")

theta_std <- theta_std %>%
  filter(pop >= pop_min & pty_dis_n >= n_min) %>%
  mutate(resp_rate = pty_dis_n / pop)

summary(theta_std$pty_dis_n)

summary(theta_std$pty_dis_n / theta_std$pop)

# Apply thresholds to draws

draws_mu <- as_tibble(x_std) %>%
  select(starts_with("mu"))

draws_theta <- as_tibble(x_std) %>%
  select(starts_with("theta")) %>%
  select(theta_std$office_index)

# Format agency estimates

theta_std <- theta_std %>%
  select(mission, dept, office, office_index, mean, sd, lb_95, ub_95) %>%
  rename(pty_dis = mean,
         pty_dis_sd = sd,
         pty_dis_lb_95 = lb_95,
         pty_dis_ub_95 = ub_95)


# Format mission estimates

mu_std <- mu_std %>%
  select(mission, mission_index, mean, sd, lb_95, ub_95) %>%
  rename(pty_dis = mean,
         pty_dis_sd = sd,
         pty_dis_lb_95 = lb_95,
         pty_dis_ub_95 = ub_95)

# Save draws for examining difference in estimates iteration by iteration

#### Format and save data for replication files ####

# Basic party disagreement data set

dis_data <- theta_std %>%
  mutate(dept = if_else(office == dept, NA_character_, dept),
         dept = if_else(dept == "The Federal Reserve", NA_character_, dept)) %>%
  rename(agency = office,
         agency_index = office_index)

save(dis_data, file = "replication_files/data/01_dis_data.RData")

readr::write_csv(dis_data, file = "replication_files/data/01_dis_data.csv")

# Mission-level disagreement

mis_data <- mu_std

save(mis_data, file = "replication_files/data/02_mis_data.RData")

readr::write_csv(mis_data, file = "replication_files/data/02_mis_data.RData")

# Draws

save(draws_theta, file = "replication_files/data/03_dis_draws.RData")

save(draws_mu, file = "replication_files/data/04_mis_draws.RData")
